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1.  Introduction 


Traditional  metal  plasticity  models,  formulated  in  terms  of  strain  rates  and  stresses  and 
incorporated  in  large-scale  numerical  analyses,  provide  useful  solutions  for  a  wide  range  of 
problems.  Details  of  the  material  microstructure  interactions  that  govern  the  deformation 
response  are  assumed  to  occur  at  length  scales  not  resolvable  by  the  simulations  and  are  captured 
implicitly  in  the  constitutive  relations.  For  example,  dependence  of  the  yield  strength  on  grain 
size  through  the  Hall-Petch  effect  can  be  incorporated  by  including  grain  size  in  the  constitutive 
model  without  tracking  the  details  of  dislocation  interactions  with  grain  boundaries. 

In  simulations  with  spatial  resolution  at  or  below  the  micron  level,  as  in  multiscale  modeling,  the 
length  scales  dictating  some  hardening  mechanisms  are  correlated  with  gradients  in  the  plastic 
strain  field.  The  torsion  experiments  of  Fleck  et  al.  (1994)  clearly  demonstrate  increased  strength 
with  decreasing  size  for  wires  10’s  of  microns  in  diameter.  The  size  effect  was  further  observed 
in  bending  (Stolken  and  Evans,  1998)  and  indentation  (Saha,  et  al.  2001),  and  in  many 
subsequent  studies.  The  strengthening  is  attributed  to  gradients  in  the  crystal  lattice  orientation. 
These  gradient  microstructures  both  store  energy  and  provide  resistance  to  further  dislocation 
motion  (Lee  et  al.,  1989;  Fleck  and  Hutchinson,  1997;  Gao  et  al.,  1999;  Baskaran  et  al.,  2010; 
Schouwenaars  et  al.,  2010).  Applying  traditional  crystal  plasticity  models  (Asaro,  1983;  Peirce  et 
al.,  1983)  to  investigate  the  size  scale  effect  in  a  multiscale  framework  will  not  be  successful 
because  the  models  are  formulated  in  terms  of  traditional  continuum  variables  of  strain  rate  and 
stress,  and  there  is  no  underlying  microstructure  length  scale  that  would  produce  a  size  effect. 

Numerous  studies  over  more  than  a  decade  have  investigated  ways  to  incorporate  a  length  scale 
into  continuum  crystal  plasticity  models.  Most  focus  on  microstructure  gradients,  as  it  is 
recognized  that  both  the  Hall-Petch  effect  and  the  size  scale  results  are  tied  to  gradients.  To  cite 
just  a  few  of  the  many  examples,  models  have  examined:  continuously  distributed  dislocations 
(Acharya,  2001),  dislocation  density  gradients  (Arsenlis  et  al.,  2004),  gradient  related  state 
variables  (Gurtin  et  al.,  2007;  Gerken  and  Dawson,  2008),  and  micro-polar  theories  (Mayeur  et 
al.,  201 1).  Many  formulations  introduce  additional  variables  into  the  solution,  which  must  be 
determined  in  concert  with  the  standard  degrees  of  freedom. 

When  using  the  traditional  crystal  plasticity  model  in  a  finite  element  code,  the  deformation  of 
neighboring  elements  is  only  connected  through  shared  nodal  displacements  and  force 
equilibrium  at  the  nodes.  Even  though  dislocations  associated  with  slip  travel  through  the 
material,  the  transmission  of  dislocations  from  one  element  to  another  is  not  represented  in  the 
model,  and  the  impact  of  coordinated  slip  on  the  deformation  field  is  not  captured.  The  indirect 
result  is  that  all  finite  element  boundaries  are  infinite  sources  and  sinks  for  dislocations. 
Advanced  crystal  plasticity  formulations,  such  as  those  based  on  lattice  orientation  gradients, 
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often  include  the  continuity  of  dislocation  flux  across  finite  element  boundaries  as  a  byproduct 
because  lattice  orientation  gradients  are  created  by  an  accumulation  of  like-signed  dislocations. 

The  goal  of  this  work  is  to  determine  if  a  simpler  approach  can  be  used  to  address  the  length 
scale  dependence  and  other  deficiencies  of  the  classical  crystal  plasticity  model.  The  starting 
point  is  simply  enforcing  dislocation  flux  continuity  across  finite  element  boundaries.  Different 
slip  conservation  approaches  were  examined  by  Hirschberger  et  al.  (201 1).  A  choice  is  made  to 
implement  the  model  in  an  explicit  finite  element  code  using  an  operator-split  approach,  where  it 
is  assumed  that  the  time  steps  are  sufficiently  small  that  the  coupled  mechanics  can  be  applied 
consecutively  over  a  time  step  rather  than  concurrently.  In  this  approach,  the  slip  in  the  elements 
can  be  treated  as  a  normal  state  variable,  whereas  it  would  be  most  effective  to  treat  it  as  a  degree 
of  freedom  in  fully  coupled,  implicit  solutions.  The  code  uses  existing  solution  variables,  and  the 
boundary  conditions  are  conceptually  straightforward.  Dislocation  flux  is  unconstrained  at  free 
surfaces,  zero  at  rigid  boundaries,  and  intermediate  at  grain  boundaries  that  are  sources  and  sinks 
for  dislocations.  The  flux  gradients  can  be  used  to  infer  the  evolution  of  dislocation  gradients  that 
lead  to  lattice  orientation  gradients. 

Initial  work  to  explore  the  feasibility  and  potential  impact  of  enforcing  dislocation  flux  between 
elements  in  a  crystal  plasticity  model  used  a  special  purpose,  explicit  dynamic  code  that 
calculated  and  enforced  an  average  slip  rate  on  element  boundaries  in  a  multi-step  constitutive 
evaluation  (Becker,  2011).  The  results  showed  the  desired  effects  of  smoothing  the  deformation 
field  and  increasing  the  hardening  rate  with  smaller  sample  sizes.  However,  the  solution 
technique  permitted  spatially  oscillatory  fields,  particularly  near  the  boundaries,  so  an  alternative 
method  for  enforcing  slip  continuity  between  elements  was  indicated. 

The  report  begins  with  an  outline  of  the  nonlocal,  gradient  crystal  plasticity  model  and  its 
implementation,  and  results  are  presented  for  single  crystal  and  polycrystal  simulations. 
Consideration  of  details  of  the  microstructure  at  micron  size  scales  prompted  exploration  and 
evaluation  of  the  semi-discrete  approach.  The  discussion  expounds  on  concerns  about  the 
adequacy  of  continuum  crystal  plasticity  models  to  represent  microstructures  at  sub-micron 
resolution. 


2.  Continuum  Model  for  Slip  Continuity 


The  goal  is  to  develop  an  enhancement  to  traditional  continuum  crystal  plasticity  models  that  can 
be  readily  incorporated  into  existing  finite  element  implementations  and  run  with  minimal 
additional  computational  overhead.  This  objective  is  facilitated  by  targeting  explicit  dynamic 
solutions  where,  because  of  the  relatively  small  time  steps,  it  is  often  possible  to  impose 
additional  constraints  through  an  operator-split  approach  in  which  the  physics  are  applied 


2 


sequentially  within  a  time  step  rather  than  concurrently.  In  the  current  context,  the  operator  split 
results  in  the  constraint  application  lagging  one  time  step  behind  its  evaluation. 

2.1  Crystal  Kinematics 

The  crystal  plasticity  model  and  implementation  used  in  the  explicit  dynamic  finite  element 
simulations  is  described  in  Becker  (2004).  It  follows  from  the  widely  used  kinematic  framework 
given  by  Asaro  (1983)  with  modifications  to  incorporate  an  equation  of  state  for  high  pressure 
applications.  The  deformation  gradient,  F,  is  notionally  decomposed  into  elastic  and  plastic 
parts,  Fe  and  Fp,  respectively; 

F  =  Fe  ■  Fp  .  (1) 

The  elastic  part  accounts  for  distortion  and  rotation  of  the  crystal  lattice,  and  the  plastic  part 
represents  slip  on  predefined  crystal  planes  and  directions  that  moves  material  but  does  not  alter 
the  underlying  crystal  lattice. 

The  velocity  gradient  obtained  from  this  kinematic  description  comprises  an  additive 
decomposition  of  an  elastic  part  and  an  inelastic  part  associated  slip  on  predefined  slip  systems: 

dx 
dx 

a=l 

s“and  ma  are,  respectively,  vectors  of  the  current  slip  directions  and  slip  plane  normals 
associated  with  each  of  the  slip  systems  (superscript  a),  and  y“are  the  corresponding  slip  rates. 
The  slip  direction  and  normal  vectors  rotate  and  distort  with  the  crystal  lattice  as  the  material 
deforms. 

2.2  Nonlocal  Model 

In  conventional  crystal  plasticity  analyses  using  an  explicit  dynamic  finite  element  code  with 
uniform  strain  elements,  the  slip  rates  are  evaluated  at  each  time  step  using  the  velocity  gradient, 
the  crystal  strength,  and  any  history  variables  associated  with  the  element.  There  is  no  direct 
connection  among  the  slip  rates  in  neighboring  elements.  In  order  to  enforce  slip  continuity 
between  elements,  a  penalty  approach  is  proposed.  Differences  in  accumulated  slip  between 
adjacent  elements  will  either  increase  or  decrease  the  resistance  to  further  slip  on  that  slip 
system.  This  requires  information  from  neighboring  elements,  making  it  a  nonlocal  method. 

The  basis  for  the  model  is  that  dislocations  move  along  the  slip  direction  from  one  element  to 
another.  The  flux  of  dislocations  crossing  an  element  boundary  for  each  slip  system  is 

0  =  Pdis  •  nf  •  (3) 

where  p%is  is  the  dislocation  density,  va  is  the  dislocation  velocity,  and  rif  is  the  outward  normal 
to  the  given  element  face.  Continuity  is  enforced  by  requiring  that  the  dislocation  flux  exiting 
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through  an  element  face  equals  the  flux  entering  the  adjacent  element  through  the  common  face. 
The  dot  product  in  equation  3  accounts  for  the  orientation  of  the  slip  system  with  respect  to  the 
element  face.  It  is  zero  if  the  slip  direction  is  parallel  to  the  face.  This  provides  the  opportunity 
for  sharp  jumps  in  slip  rates  across  parallel  slip  planes  while  enforcing  continuity  along  slip 
planes. 

The  dislocation  density,  velocity,  and  the  Burgers  vector,  b,  are  related  to  the  continuum  slip  rate 
by  Orowan’s  equation 

Ya  =  Pais  va  b.  (4) 

For  a  shared  element  face,  and  assuming  that  the  Burgers  vector  is  constant  and  that  the  slip 
directions  are  closely  aligned  across  the  interface,  the  continuity  error  in  the  accumulated 
dislocation  flux  between  elements  can  be  approximated  as 

[(.Ya)neighbor- (,Ya)e\sa  -nf  =  Ef  (5) 

The  subscripts  e  and  neighbor  on  the  accumulate  slip,  ya,  denote,  respectively,  the  element  of 
interest  and  the  neighbor  sharing  the  face.  The  approach  taken  is  to  penalize  the  error.  It  is  not 
necessary  to  assume  that  the  Burgers  vectors  are  equal  or  that  the  slip  directions  in  the  adjacent 
elements  are  aligned,  but  these  simplifications  are  made  for  expedience  and  coding  clarity  in  this 
initial  implementation. 

Application  of  equation  5  in  a  penalty  method  requires  a  cumulative  result  over  all  element  faces. 
This  is  complicated  by  the  changing  algebraic  sign  of  the  dot  product  for  different  element  faces 
and  the  50%  probability  that  the  accumulated  slips  are  negative.  The  sign  of  the  dot  product  is 
irrelevant  to  the  physical  problem,  so  the  absolute  value  is  taken.  The  algebraic  sign  of  the  slip  is 
also  unimportant,  and  the  magnitude  is  obtained  by  multiplying  by  the  sign  of  the  slip.  This  is 
preferred  to  using  the  absolute  values  of  the  slips  in  cases  where  the  sign  of  the  slip  is  different  in 
the  neighboring  element.  With  these  adjustments,  the  driving  force  for  the  penalty  method, 
equation  5,  is  summed  as 

N  faces 

- T  y  7~SiSKYa)[(ya)neighbor- (Ya)e]\sa  -nf\  =  Ea  (6) 

i  r\  f=1  j 

Lf  is  the  distance  between  element  centroids  associated  with  the  particular  element  face.  It  is 
intended  to  provide  a  larger  penalty  if  the  slip  difference  occurs  on  a  smaller  spatial  grid.  Ea  in 
equation  6  is  a  weighted  sum  of  unsigned  slip  gradients  on  a  slip  system  for  a  given  element.  If 
the  value  is  positive,  slip  is  deficient  in  the  element  and  additional  slip  is  promoted.  If  the  value 
is  negative,  further  slip  is  impeded.  It  is  significant  that  a  constant  gradient  results  in  face 
contributions  that  sum  to  zero.  Hence,  a  constant  gradient  is  not  suppressed.  The  excess 
dislocations  associated  with  the  gradient  are  assumed  to  be  uniformly  distributed. 
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A  closer  look  at  equation  6  for  the  special  case  of  a  slip  direction  aligned  with  a  regular 
rectangular  mesh  reveals  a  serendipitous  connection  to  published  gradient  theories.  Assuming 
that  the  signs  of  the  slips  are  the  same  and  performing  the  summation. 


Subscripts  e—1  and  e+1  denote  elements  on  either  side  of  element  e.  Dividing  equation  7  by  Lf, 
this  becomes  the  finite  difference  form  for  half  of  the  second  derivative  of  the  slip  along  the  slip 
direction.  The  second  derivative  appears  in  the  micro  force  balance  in  several  gradient 
formulations  (e.g.,  Bittencourt  et  al.,  2003;  Gurtin  et  al.,  2007).  Hence,  with  a  small  coding 
change,  the  implementation  can  approximate  a  simplified  version  of  an  established  gradient 
model.  Because  of  this  connection  to  established  models,  and  since  it  gives  a  nonlinear  length 
scale  dependence,  the  modified  relation 


is  used  for  the  nonlocal  simulations. 


2.3  Crystal  Flow  Strength 

The  rate  dependent  crystal  constitutive  model  relates  the  loading  on  each  slip  system  to  the  slip 
rate.  The  driving  force  for  slip,  Ta,  is  the  resolved  shear  stress  on  the  slip  system. 


where  a  is  the  Cauchy  stress.  A  simple  power  law  rate  model  is  used  here  to  determine  the  slip 
rate  for  the  applied  loading. 

/  Ta  \1,m 

7“  =  To  (  ^  j  sign(ra )  .  (10) 


A  low  rate  exponent,  m  =  0.005,  is  chosen  to  provide  nearly  rate-independent  behavior  while 
still  smoothing  the  transition  from  elastic  to  plastic  response  at  the  slip  system  flow  strength,  ga. 

The  flow  strength  is  adjusted,  based  on  the  nonlocal  contribution,  to  increase  or  reduce  the 
resistance  to  continued  deformation 

ga  =  g  o  —  gpi  tanh  [b2Ea  —  ^  (ii) 


Here,  p1  and  p2  are  dimensionless  parameters,  p,  is  some  average  shear  modulus  for  the  crystal, 
and  b  is  the  Burgers  vector.  The  base  crystal  strength,  g$,  could  be  a  function  of  slip  to  capture 
strain  hardening,  but  it  is  assumed  constant  in  these  analyses  to  simplify  interpretation  of  the 
results.  The  hyperbolic  tangent  function  is  used  to  cap  the  influence  of  the  nonlocal  term.  The 
contribution  is  approximately  linear,  ~pp2b2Fa,  while  the  argument  is  small,  and  it  is  capped  at 
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a  constant  value  of  ±pp1.  The  penalty  parameters  are  chosen  as  px  =  5  x  10  4  and  p2  = 

10000.  In  all  of  the  nonlocal  continuum  simulations:  the  base  flow  strength  is  =  33.33  MPa; 

Q 

the  density  is  p0  =2.7  g/cm  ;  the  Burgers  vector  is  b  =  0.3  nm;  and  shear  modulus  is  p  = 
25,650  MPa. 

2.4  Crystal  Geometry  and  Boundary  Conditions 

An  idealized  two-dimensional  crystal  geometry  is  used  for  these  analyses.  The  crystal  consists  of 
three  slip  systems  set  in  an  equilateral  triangle  configuration.  This  allows  a  multitude  of  slip 
modes,  but  redundant  slip  is  unlikely.  In  the  rate  independent  limit,  it  is  not  possible  to  have  slip 
on  all  three  systems  simultaneously,  and  the  rate  sensitivity  is  too  low  to  enable  redundant  slip 
activation  in  these  simulations. 

Constant  strain  quadrilateral  elements  with  hourglass  control  (Flanagan  and  Belytschko,  1981) 
are  used  for  all  of  the  simulations.  Elements  in  which  all  four  faces  are  adjacent  to  an  element  of 
the  same  initial  orientation  are  interior  elements,  and  gradients  are  computed  directly  as  indicated 
in  equation  8.  Elements  with  fewer  than  four  faces  contacting  regions  of  the  same  orientation  are 
either  on  grain  boundaries  or  model  boundaries.  The  face  is  flagged  for  these  elements,  and  a 
parameter  is  checked  to  see  whether  it  is  treated  as  a  zero  flux  boundary  or  a  free  boundary  with 
no  slip  impedance.  For  the  non-interior  elements  with  zero-flux  boundaries,  ghost  elements  with 
opposite  slip  are  assumed  across  the  flagged  faces.  Grain  boundaries  and  surfaces  with  applied 
boundary  conditions  are  treated  in  this  manner.  For  free  surfaces  and  periodic  boundaries,  the 
ghost  element  across  the  flagged  face  is  set  with  the  same  slip  so  that  these  interfaces  do  not 
contribute  to  the  gradient. 

2.5  Finite  Element  Implementation 

The  model  was  implemented  in  the  large-scale  parallel,  explicit  finite  element  code  ALE3D 
(2012).  The  crystal  plasticity  constitutive  model  existed  previously  (Becker,  2004)  and  the 
strengthening  terms  in  equation  1 1  due  to  the  nonlocal  effects  were  a  straightforward  addition. 

The  nonlocal  computations  occur  outside  of  the  material  model  when  all  history  variables  are  at 
a  consistent  state.  Since  it  is  an  Arbitrary  Lagrange-Eulerian  code  that  moves  material  through 
the  mesh,  many  features  and  functions  are  already  in  place  for  the  nonlocal  calculations.  Lists  of 
neighboring  elements  and  shared  faces  exist,  as  do  functions  accessing  history  variables  in  the 
adjacent  elements.  For  parallel  computations,  the  problem  is  subdivided  into  domains  that  reside 
on  separate  processors.  Elements  that  are  at  the  boundaries  of  these  domains  have  neighboring 
elements  that  reside  on  other  processors.  The  information  from  the  neighboring  elements 
residing  on  other  processors  is  carried  locally  in  ghost  elements,  so  all  of  the  data  needed  to 
perform  the  nonlocal  calculations  of  equation  8  are  available  locally  on  each  processor.  A 
communication  call  prior  to  the  nonlocal  calculations  assures  that  all  of  the  information  is 
current.  A  significant  amount  of  bookkeeping  and  data  juggling  is  required,  but  the  computations 
in  equation  8  are  straightforward. 
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3.  Continuum  Finite  Element  Results 


The  effect  of  the  slip  continuity  constraint  (slip  gradient)  is  evaluated  on  two  configurations, 
each  at  multiple  length  scales.  All  simulations  are  two-dimensional.  The  first  configuration  is 
simple  shear  of  a  single  crystal  with  one  of  the  slip  planes  initially  aligned  orthogonal  to  the 
shear  direction.  This  creates  single  slip  conditions  for  a  straightforward  evaluation  of  the  model. 
The  second  configuration  is  a  polycrystal  constructed  from  regular  hexagons.  The  orientation  of 
the  crystal  lattice  for  each  grain  is  random. 

3.1  Single  Crystal 

Single  crystal  calculations  were  run  at  four  size  scales  using  a  20  x  100  mesh  of  square  elements 
(figure  1).  Velocity  boundary  conditions  are  applied  to  the  upper  and  lower  surfaces  to  shear  the 
top  of  the  crystal  to  the  right.  Initial  velocities  of  all  interior  nodes  are  prescribed  consistent  with 
simple  shear  to  eliminate  ringing  as  the  explicit  dynamic  calculation  starts.  Periodic  boundary 
conditions  are  applied  to  the  lateral  surfaces  to  mimic  an  infinitely  wide  crystal.  The  heights  of 
the  single  crystals  simulated  were  50,  5,  1,  and  0.5  pm,  and  the  width  of  the  simulation  box  was 
20%  of  the  height  in  each  case.  Although  the  width  is  irrelevant  with  the  periodic  boundary 
conditions,  multiple  elements  are  used  across  the  width  to  demonstrate  that  the  boundary 
conditions  are  applied  properly.  Slip  transmission  is  restrained  on  the  upper  and  lower 
boundaries,  and  slip  transmission  is  unimpeded  on  the  lateral,  periodic  boundaries. 


Velocity  boundary  sw 
condition 

Lattice 

Periodic  boundary 
conditions  \H| 

/  orientation 

<r 

Figure  1 .  Initial  configuration  for  the  single  crystal  and  the  crystal  lattice  orientation. 
The  bottom  is  fixed  and  the  top  is  moved  to  the  right.  Periodic  boundary 
conditions  are  applied  coupling  the  left  and  right  hand  sides. 


3.1.1  Nonlocal  Strength 

The  nonlocal  strength  contribution  to  equation  1 1  is  shown  in  figure  2  for  the  four  crystal  sizes 
and  at  two  shear  strains.  The  largest  crystal  is  on  the  left  and  the  smallest  is  on  the  right.  The  top 
row  shows  the  distribution  at  a  shear  strain  of  0.03  and  bottom  shows  the  results  at  an  average 
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shear  strain  of  0.05.  The  color  scales  are  different  for  each  crystal  size,  but  the  scales  at  each 
crystal  size  are  the  same  at  0.03  and  0.05  shear  strains  to  highlight  the  evolution. 


Figure  2.  Distribution  of  the  nonlocal  contribution  to  the  strength  on  the  slip  system  aligned  vertically  in  the 
crystal,  for  the  crystal  thicknesses  indicated,  at  shear  strains  of  0.03  and  0.05. 

The  magnitude  of  the  nonlocal  strength  contribution  and  the  relative  depth  that  the  distribution 
penetrates  from  the  surfaces  increases  as  the  crystal  thickness  is  decreased.  The  sharpest 
gradients  are  expected  at  the  crystal  surfaces  where  the  slip  transmission  is  impeded.  For  the 
50-pm  crystal  the  strength  is  increased  only  within  a  few  elements  of  the  surface,  corresponding 
to  a  few  microns.  The  central  portion  of  the  crystal  sees  no  gradient  or  strengthening  effect,  even 
as  the  strain  increases  from  0.03  to  0.05.  The  boundary  layer  also  appears  to  penetrate  a  few 
microns  in  the  5-pm-thick  crystal  simulation.  However,  a  smaller  portion  of  the  crystal  is  nearly 
free  of  gradient  effects  for  this  smaller  crystal.  At  yet  smaller  crystal  thicknesses,  the  gradient 
effect  penetrates  the  full  crystal  thickness,  and  the  evolution  with  increasing  deformation  is 
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evident.  The  strength  at  the  center  is  elevated  significantly  by  a  strain  of  0.03,  and  it  continues  to 
increase  with  further  deformation. 

The  minimum  and  maximum  values  are  indicated  in  each  of  the  plots.  Except  for  the  smallest 
crystal,  the  difference  between  the  maximum  and  minimum  increases  as  the  crystal  size 
decreases  and  the  difference  also  increases  with  increasing  strain.  The  increase  with  strain 
indicates  that  the  gradient  is  still  evolving  at  a  shear  strain  of  0.05.  The  trends  are  different  for 
the  0.5 -pm  crystal  because  the  gradient  strengthening  is  beginning  to  run  up  against  the  cap  set 
by  the  hyperbolic  tangent  function  in  equation  11.  With  continued  deformation,  the  gradient 
strengthening  is  becoming  more  uniform,  albeit  at  a  higher  level. 

3.1.2  Stress  State 

The  momentum  balance  in  the  vertical  direction  (y-direction)  dictates  that  the  y-gradient  in  the  y- 
direction  stress  component  is  balanced  by  the  horizontal  (x-direction)  gradient  in  the  shear  stress. 
Since  the  periodic  boundary  conditions  require  that  all  horizontal  gradients  are  zero,  the  stress  in 
the  y-direction  should  be  constant.  The  magnitude  of  the  y-direction  stress  is  not  determined  by 
the  momentum  equation,  just  that  it  is  constant.  The  calculations  show  a  constant  y-direction 
stress  to  five  significant  digits. 

More  important  for  current  purposes  is  the  momentum  balance  in  the  horizontal  direction.  The 
lack  of  stress  gradients  in  the  horizontal  direction  requires  the  shear  stress  to  be  constant  through 
the  thickness.  The  simulations  show  that  the  shear  stress  is  constant  to  five  significant  digits. 
While  the  x-direction  stress  must  be  constant  in  the  x-direction,  the  symmetry  conditions  and 
momentum  equations  provide  no  further  constraints  restricting  its  gradient  in  the  y-direction. 

3.1.3  Slip  Rate 

The  normalized  slip  rates  corresponding  to  the  configurations  in  figure  2  are  shown  in  figure  3. 
Again,  the  largest  crystal  is  on  the  left  and  the  smallest  is  on  the  right.  The  top  row  contains 
results  at  a  shear  strain  of  0.03,  and  bottom  row  shows  the  normalized  slip  rate  at  a  shear  strain 
of  0.05.  The  plots  are  normalized  by  the  applied  shear  strain  rate  so  that  a  value  of  1.0  would 
indicate  a  uniform  shear.  Since  the  shear  stress  and  the  reference  strength  in  equation  1 1  are  both 
constant,  the  slip  rate  is  approximately  related  to  the  gradient  term  through  equations  10  and  11. 
Second-order  factors  influencing  the  slip  rate  include  the  change  in  lattice  orientation  and  non¬ 
zero  components  of  the  x-direction  and  y-direction  stresses  that  modify  the  resolved  shear  stress. 
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Figure  3.  Distribution  of  the  normalized  slip  rate  for  single  crystals  of  the  indicated  thickness  and  at  shear  strains 
of  0.03  and  0.05.  The  slip  rates  are  normalized  by  the  applied  shear  rate. 

The  slip  rate  follows  the  same  trends  as  the  nonlocal  hardening  contribution.  The  variation  in  slip 
rate  is  greater  for  the  smaller  crystals;  and,  except  for  the  smallest  crystal,  the  differences  are 
greater  with  increased  deformation.  As  a  result  of  the  boundary  conditions,  slip  rates  are  low  at 
top  and  bottom  boundaries  compared  to  the  center  regions.  This  results  in  the  sigmoidal 
deformed  shapes.  More  severe  differences  in  slip  rate  result  in  greater  deviation  from  a  linear 
shear  deformation  profile.  The  slip  rates  for  the  0.5-pm  crystal  become  more  uniform  at  the 
higher  deformation  because  the  gradient  term  is  capped  by  the  hyperbolic  tangent  function.  The 
strength  is  more  uniform,  which  results  in  a  more  uniform  slip  rate. 

3.1.4  Stress  Strain  Response 

The  shear  stress-shear  strain  responses  for  the  various  crystal  thicknesses  are  plotted  in  figure  4. 
The  curves  are  identical  through  the  linear  elastic  regime,  and  all  yield  at  approximately  the  same 
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stress,  approximately  34  MPa.  The  applied  shear  rate  is  50  times  the  reference  shear  rate;  and, 
accounting  for  the  strain  rate  sensitivity,  the  apparent  yield  strength  is  calculated  from  equation 
10  to  be  33.99  MPa  rather  than  the  reference  shear  strength  of  33.33  MPa. 


Shear  Strain 


Figure  4.  Shear  stress-shear  strain  response  predicted  for  the  four  crystal  thicknesses. 

The  stress  shows  the  expected  trend  of  increasing  strength  at  the  smaller  crystal  sizes.  This  is 
directly  related  to  the  nonlocal  strength  in  figure  2.  The  nonlinear  dependence  on  specimen 
dimensions  is  due  to  dividing  by  the  element  size  squared  in  equation  8.  The  strain  hardening 
rate  is  fairly  consistent  with  increasing  strain  for  the  three  larger  crystal  sizes  but  not  for  the 
smallest.  Stress  in  the  0.5-pm  crystal  peaks  as  the  hyperbolic  tangent  function  places  a  cap  on  the 
nonlocal  strength  contribution.  This  is  also  consistent  with  the  results  presented  in  figures  2  and 
3.  A  final  observation  from  figure  4  is  the  kink  that  is  most  evident  in  the  larger  two  specimens 
near  a  strain  of  0.02.  This  marks  the  transition  from  single  slip  at  lower  strains  to  slip  on  two  slip 
systems  at  larger  strains.  As  the  crystal  lattice  rotates  and  stresses  build  in  the  x  and  y  directions, 
the  crystals  are  able  to  accommodate  the  deformation  more  easily  with  multiple  active  systems. 
Due  to  the  angle  of  the  slip  plane,  slip  constraints  at  the  boundaries  are  not  as  severe  for  the 
second  slip  system,  so  the  strain  hardening  rate  is  reduced. 

3.1.5  Mesh  Refinement 

The  effect  of  mesh  resolution  on  the  solution  is  investigated  by  rerunning  the  5-pm-thick 
simulation  using  twice  as  many  elements  in  each  direction.  The  results  from  the  40x200  mesh  are 
shown  along  side  of  the  20x100  mesh  results  in  figure  5.  Other  than  the  expected  differences  in 
smoothness  of  the  fields,  the  nonlocal  stress  and  slip  rate  distributions  do  not  appear  to  be 
influenced  significantly  by  halving  the  mesh  size. 
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Figure  5.  Comparison  of  nonlocal  stress  and  the  slip  rate  for  the  20x100  and  40x200  element  simulations  of 
the  5-pm-thick  single  crystal.  The  results  are  shown  at  a  0.05  shear  strain. 

For  a  more  quantitative  assessment,  the  applied  shear  stress  computed  at  the  crystal  boundaries  is 
2.3%  lower  for  the  finer  mesh.  Part  of  the  difference  may  be  due  to  the  lower  energy  solution 
expected  with  an  increased  number  of  degrees  of  freedom,  and  the  remainder  can  be  attributed  to 
discretization  error  associated  with  the  nonlocal  computations  and  the  operator  split  algorithm. 
The  time  step  was  also  a  factor  of  two  lower  in  the  fine  mesh  calculation  due  to  the  dependence 
of  the  Courant  stable  time  step  on  the  mesh  size.  The  smaller  time  step  improves  the  accuracy  of 
the  operator  split  integration. 

3.1.6  Time  Step  Instability 

The  simulation  of  the  0.5-pm-thick  crystal  experienced  numerical  instabilities  when  the  strain 
increment  per  step  was  too  great.  This  is  thought  to  be  associated  with  the  operator  split  where 
the  strength  increase  from  the  slip  gradients  creates  a  driving  force  that  is  too  large  and  over¬ 
corrects  the  slip  rate.  Specifically,  with  a  strain  increment  of  1.54x10  7  per  time  step,  the  slip  rate 
for  the  0.5-pm  crystal  was  erratic  and  non-zero  only  in  scattered,  isolated  elements.  These 
isolated  regions  of  slip  occurred  briefly  and  died  out  quickly  as  deformation  proceeds,  and 
eventually  strain  was  incremented  in  the  entire  domain,  albeit  unevenly.  When  the  time  step  was 
reduced  by  a  factor  of  two,  such  that  the  strain  increment  per  step  was  7.43x1 0-8,  the  calculation 
was  well  behaved.  The  results  in  figures  2  through  4  were  run  with  a  strain  increment  of 
3.853x10  8  to  be  certain  that  the  time  step  was  small  enough  to  suppress  the  instability.  The 
calculations  for  the  larger  crystal  sizes  experienced  less  gradient  hardening,  and  they  were  run  at 
the  Courant  stable  time  step  without  any  additional  time  step  controls. 
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3.2  Polycrystal 


An  idealized  polycrystal  was  created  by  filling  a  rectangular  region  with  regular  hexagonal 
grains  (figure  6a).  Simple  shear  boundary  conditions  were  applied  by  prescribing  velocities  to 
the  nodes  on  the  upper  and  lower  surfaces.  Periodic  boundary  conditions  were  applied  on  the  left 
and  right  surfaces.  The  orientation  of  the  triangular  crystal  lattice  in  each  of  the  grains  was 
random,  and  the  rotation  angle  from  the  reference  orientation  is  indicated  on  the  plot.  The  half 
grains  at  the  same  height  on  the  left  and  right  of  the  model  region  were  given  the  same 
orientation  to  facilitate  application  of  periodic  boundary  conditions.  As  with  the  single  crystal 
simulations,  the  initial  velocity  of  all  interior  nodes  was  prescribed  to  eliminate  ringing  from 
abrupt  imposition  of  boundary  conditions. 


The  default  inter-element  slip  rate  condition  for  all  elements  is  that  any  element  face  touching 
another  grain  will  have  restricted  slip.  This  is  imposed  by  assuming  that  a  ghost  element  across 
the  interface  has  equal  and  opposite  slip  in  equation  8.  The  restricted  slip  condition  is  enforced 
on  the  upper  and  lower  surfaces  and  on  grain  boundaries,  including  those  grain  boundaries  on  the 
periodic  surfaces.  The  half  crystals  on  the  periodic  surfaces  are  treated  differently;  the  element 
across  the  interfaces  is  assumed  to  have  the  same  slip.  This  is  not  a  truly  periodic  condition,  but  a 
data  structure  identifying  periodic  neighboring  elements  is  not  yet  available. 


Three  model  sizes  are  investigated:  160\/3  /mi  x  300  /mi,  l6^f3~|m\  x  30  /mi,  and 
1.6a/3  /an  x  3  /t m.  All  use  the  same  mesh  configuration,  scaled  to  give  the  appropriate 
dimensions.  Each  of  the  88  hexagons  was  discretized  by  an  identical  mesh  of  2112  quadrilateral 
elements  (figure  6b).  Nodes  are  shared  along  the  grain  boundaries,  so  the  deformation  is 
continuous  throughout. 


Figure  6.  Grain  structure  (a)  and  finite  element  mesh  and  (b)  for  the  polycrystal  simulations. 
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3.2.1  Slip  Rate 


Slip  rates  normalized  by  the  applied  shear  rate  are  shown  in  figure  7  for  the  three  polycrystal 
sizes.  At  the  larger  crystal  size,  the  strain  rate  localizes  into  well-defined  bands.  The  majority  of 
the  deformation  is  carried  by  two  horizontal  bands  with  some  scattered  activity  in  the  central 
region  of  the  model.  Blocks  of  grains  appear  to  remain  elastic  while  localized  shear  along  the 
horizontal  and  vertical  bands  accommodates  the  deformation  between  neighboring  blocks.  The 
bands  do  not  follow  the  grain  boundaries,  but  many  are  associated  with  grain  boundary  triple 
points. 


At  the  intermediate  polycrystal  size,  the  nonlocal  slip  constraint  diffuses  the  deformation.  The 
slip  bands  are  still  fairly  well  defined,  but  the  peak  strain  rates  are  not  as  high  and  regions  of 
nearly  elastic  behavior  are  smaller  and  less  well  defined.  The  slip  rates  in  the  3.3-pm-thick 
polycrystal  are  considerably  more  diffuse  and  the  material  near  the  highly  constrained  top  and 
bottom  boundaries  has  the  lowest  strain  rates.  Grain  outlines  are  evident  as  the  slip  rate  tends  to 
be  high  or  low  at  the  grain  boundaries,  and  the  color  contrast  across  the  boundaries  accentuates 
them. 


Figure  7.  Normalized  slip  rate  distribution  for  simple  shear  deformation  of  idealized  polycrystals  with  heights  of 
330,  33,  and  3.0  pm.  The  color  scale  is  the  same  for  all  three  plots. 

3.2.2  Nonlocal  Strengthening 

The  nonlocal  strengthening  associated  with  the  gradients  is  shown  in  figure  8  for  all  three  slip 
systems  and  the  three  crystal  sizes.  The  color  scale  in  each  row  is  the  same  so  that  the  magnitude 
of  the  effect  of  the  slip  systems  can  be  compared.  The  scales  are  different  for  each  crystal  size  as 
the  strengthening  is  much  greater  in  the  smaller  model  region.  The  color  scale  on  for  the  330-mm 
polycrystal  is  set  to  a  relatively  low  value  of  0.5  MPa,  and  even  then,  the  gradient  contribution  is 
only  evident  at  the  grain  boundaries  or  near  the  most  highly  shear  regions  shown  in  figure  7.  The 
impact  on  the  solution  is  minor. 
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Figure  8.  Nonlocal  strength  contributions  on  the  three  slip  systems  for  polycrystal  model  sizes  of  330,  33,  and 
3.3  pm.  The  color  scales  are  consistent  within  each  row. 

The  nonlocal  strength  distribution  in  the  3 3 -(am  polycrystal  appears  in  the  grain  interiors  as  well 
as  at  grain  boundaries.  The  strong  interior  features  are  associated  with  stress  concentrations  at 
grain  boundary  triple  points,  and  most  correspond  to  elevated  slip  activity  in  figure  7.  Many  of 
the  grain  boundaries  show  strengthening  on  one  side  and  softening  on  the  other.  These 
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correspond  to  increasing  slip  when  approaching  grain  boundaries  for  strengthening  and 
decreasing  slip  when  approaching  the  boundaries  for  softening. 

The  nonlocal  strengthening  is  quite  prominent  in  the  3.3-pm  polycrystal.  As  with  the 
intermediate  size  polycrystal,  the  sign  of  the  gradient  effect  is  often  flipped  across  the  grain 
boundaries.  For  most  grains,  the  gradient  is  strongest  at  the  grain  boundaries  and  decays  toward 
the  grain  center.  However,  there  are  a  few  notable  grains  where  the  peak  values  are  in  the 
interiors.  These  correspond  to  locations  of  intersecting  slip  activity  in  figure  7.  The  hyperbolic 
tangent  function  causes  the  gradient  effect  to  saturate  at  a  level  between  12  and  13  MPa.  It  is  also 
notable  that  the  gradient  strengthening  occurs  on  only  two  of  the  three  slip  systems.  This  reflects 
the  lack  of  redundant  slip  for  the  idealized  crystal.  Only  two  slip  systems  are  active  at  any  time. 

3.2.3  Polycrystal  Stress-strain  Behavior 

The  shear  stress-shear  strain  response  for  the  three  polycrystal  sizes  is  presented  in  figure  9.  As 
with  the  single  crystals,  the  stress  is  higher  for  the  smaller  polycrystals.  The  nonlinearity  with 
length  scale  is  also  clearly  evident.  There  are,  however,  two  important  distinctions  from  the 
single  crystal  results.  The  first  is  that  the  initial  yield  point  varies  with  crystal  size,  whereas  it  did 
not  for  the  single  crystal  simulations.  This  is  thought  to  be  related  to  the  single  crystals  yielding 
throughout  simultaneously,  while  the  polycrystal  yields  gradually  and  may  build  up  local  slip 
gradients  before  the  macroscopic  yield  is  evident. 


Shear  Strain 


Figure  9.  Shear  stress  strain  response  for  three  different  size  scales 
of  idealized  polycrystals. 

The  other  notable  difference  is  that  the  curves  are  not  smooth.  This  could  result  from  a 
combination  of  the  evolution  of  the  crystal  lattice  orientation,  evolution  of  the  slip  gradients,  and 
wave  propagation  in  the  explicit  dynamic  calculation.  The  change  in  lattice  orientation  is  shown 
in  figure  10  for  the  largest  and  smallest  size  scales.  In  the  330-pm  polycrystal,  where  the  strain 
localization  is  more  pronounced,  lattice  reorientation  is  also  localized.  The  local  geometric 
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softening  facilitates  the  shear.  In  contrast,  due  to  the  slip  continuity  and  gradient  constraints,  the 
lattice  reorientation  in  the  3.3-pm  polycrystal  is  smoothed  over  a  larger  region  relative  to  the 
grain  size,  and  the  lattice  within  the  grains  rotates  nearly  uniformly.  A  larger  portion  of  the 
polycrystal  had  to  realign,  which  takes  longer  and  results  in  a  greater  load  excursion  before  it 
settles  into  a  nearly  steady  shear  mode. 


Figure  10.  Change  in  crystal  lattice  orientation,  in  degrees,  at  0.025  shear  strain  for  the  330-  and  3.3-pm  high 
polycrystals. 


4.  Considerations  for  Discreteness  of  Dislocations 


The  motivation  for  imposing  dislocation  flux  constraints  is  that  dislocations  are  discrete  entities 
that  propagate  from  one  element  to  the  next  as  part  of  the  slip  process.  Another  aspect  of  the 
dislocation  discreteness  is  their  spacing,  which  is  typically  quantified  in  terms  of  the  dislocation 
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density.  For  well-annealed  metals,  a  typical  dislocation  density  is  10  -  10  cm  ;  at  a  few 
percent  deformation,  this  increases  to  10s  -  109  cm-2;  and  for  a  very  heavily  deformed 
polycrystal,  the  dislocation  density  is  in  the  neighborhood  of  1011  cm  2  (Hull  and  Bacon  1984). 
Table  1  lists  the  element  areas  for  each  of  the  four  simulations  in  section  3.1  and  the  higher 
values  of  dislocation  density  for  well-annealed,  lightly  deformed,  and  heavily  deformed 
polycrystalline  metals.  From  these  values,  an  average  number  of  dislocations  enclosed  by  an 
element  is  calculated. 
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Table  1.  Average  number  of  dislocations  per  element  expected  in  the  simulations  of  section  3.1  for  well -annealed, 
lightly  deformed,  and  heavily  deformed  metals. 


50  jLim  Crystal 
Elength  =  0.5  rim 
Earea  =  0.25  Mm2 

5  jum  Crystal 
^length  =  0.05  JLim 
Earea  =  0.0025  jLim2 

1  jum  Crystal 

Eie„gth  =  0.01  ym 
Earea  =  0.0001  (UTl2 

0.5  jum  Crystal 
Eiength  =  0.005  Mm 
Earea  =  0.000025  Jill/ 

107  cm  2  (10  1  |um  2) 

1/40 

1/4000 

10“5 

2.5x10^ 

109  cm-2  (101  jam  2) 

2.5 

1/40 

1/1000 

1/4000 

to11  cm"2  (to3  |jm-2) 

250 

2.5 

1/10 

1/40 

The  entries  in  table  1  that  are  less  than  1.0  indicate  that  not  every  element  will  contain  a 
dislocation.  For  example,  1/4000  means  that  one  out  of  every  4000  elements  can  be  expected  to 
contain  a  dislocation.  An  implicit  assumption  in  continuum  crystal  plasticity  models  is  that  the 
dislocation  content  in  the  elements  is  sufficient  for  slip  to  be  smooth  and  continuous.  It  is  clear 
that  these  conditions  are  not  met  for  a  well-annealed  metal  using  any  of  the  meshes  in  section 
3.1,  since  not  every  element  would  contain  even  one  dislocation.  Using  the  discretization 
provided  by  the  50-pm-thick  crystal  simulation,  a  sufficient  number  of  dislocations  would  be 
represented  within  each  element  when  the  crystal  is  heavily  deformed,  but  not  in  the  deformation 
leading  up  to  that  state.  Considering  that  dislocations  are  usually  not  uniformly  distributed, 
element  sizes  of  a  few  microns  may  be  necessary  to  assure  a  sufficient  number  of  dislocations 
per  element  for  a  proper  continuum  crystal  model  representation. 

Several  aspects  of  crystal  deformation  cannot  be  represented  accurately  by  a  continuum  model  if 
elements  are  sparsely  populated  by  dislocations.  One  difficulty  is  that  only  elements  containing 
dislocations  can  slip,  and  other  elements  must  either  deform  elastically  or  nucleate  additional 
dislocations.  A  second  feature  not  well  represented  is  the  discreteness  of  slip.  A  dislocation  that 
passes  through  an  element  creates  a  slip  increment  Ay  =  b/Eiength.  For  a  Burgers  vector  length 

of  b  =  0.3xl0-3  pm  and  an  element  size  of  E[ength  =  0.5  pm,  the  slip  must  occur  in  increments  of 
0.6  xlO  3  rather  than  in  the  arbitrarily  small  increments  permitted  by  the  continuum 
representation.  The  discrete  strain  increment  as  a  dislocation  moves  from  one  element  to  the  next 
causes  a  commensurate  increment  in  stress.  This  stress  jump  can  be  comparable  to  the  yield 
strength.  A  final  aspect  of  crystal  plasticity,  which  cannot  be  represented  properly  if  the  element 
size  is  less  than  the  dislocation  spacing,  is  dislocation  interactions  leading  to  strengthening. 
Hardening  and  gradient  effects  are  a  result  of  interactions  among  dislocations.  If  the  dislocations 
are  sparse  in  the  mesh,  these  interactions  must  be  accounted  for  explicitly  rather  than  implicitly 
in  a  hardening  function. 

4.1  Semi-discrete  Dislocation  Model 

In  an  attempt  to  push  continuum  finite  element  simulations  to  smaller  length  scales  where 
dislocations  are  sparse  within  the  elements,  a  semi-discrete  model  was  developed.  It  is  run  within 
a  standard  explicit-dynamic  finite  element  framework  that  is  described  in  Becker  (201 1).  The 
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single  slip  constitutive  model  follows  the  traditional  formulation  with  three  major  modifications: 
(1)  only  elements  that  contain  dislocations  or  dislocation  sources  can  slip;  (2)  the  slip  increment 
is  quantized  in  terms  of  the  Burgers  vector  and  element  size;  and  (3)  elements  designated  to 
contain  dislocation  nucleation  sources  have  a  reduced  flow  strength.  In  addition  to  these  slip 
model  modifications,  the  code  tracks  dislocations  moving  from  one  element  to  another,  and  it 
also  tracks  the  total  number  of  dislocations  that  have  traversed  each  element. 


The  quantization  of  the  slip  increment  is  created  by  starting  with  the  standard  continuum  slip 
increment,  A Y continuums  and  using  it  to  create  an  integer  representing  the  number  of  Burgers 
vector  lengths  the  dislocation  traversed  during  increment. 


NBurg  =  int 


a,,  length  j 

^ Y COTltiTllilLTTl  l  ) 


(12) 


The  square  of  the  dimensionless  size  scale  loosely  accounts  for  the  step-by-step  dislocation 
progression  across  the  element.  The  slip  increment  used  in  the  stress  update  is  then  calculated  by 

A Yupdate  =  N Burg  )  (13) 

\^length  J 

Similar  integer  arithmetic  is  used  to  determine  when  a  dislocation  has  completely  traversed  an 
element  and  will  be  transferred  to  a  neighbor. 

The  crystal  is  assumed  to  generate  dislocations  dipoles  at  fixed  nucleation  sites.  These 
dislocation  nucleation  sites  are  assigned  randomly  at  the  start  of  the  simulation  to  a  small 
fraction  of  the  elements  using  a  random  number  generator.  The  yield  strength  in  these  elements 
is  also  set  to  vary  randomly  between  25%  and  75%  of  the  yield  stress.  If  the  resolved  shear  stress 
in  these  elements  exceeds  the  local  yield  strength,  the  crystal  slips  by  moving  dislocations  from  a 
dipole  in  opposite  directions  along  the  slip  plane.  Eventually  these  dislocations  move  to  the 
neighboring  elements.  At  this  time  the  dislocations  are  transferred  to  the  neighboring  elements, 
and  a  new  dislocation  dipole  initiates  at  the  nucleation  site.  Elements  with  dislocations,  but  not 
nucleation  sites,  slip  when  the  resolved  shear  stress  exceeds  the  full  50-MPa  yield  strength. 
Elements  without  dislocations  deform  elastically. 

The  semi-discrete  dislocation  model  is  assessed  through  simulations  using  square  elements  with 
sizes  1.25,  0.125,  and  0.025  pm.  The  element  size  is  important  as  it  is  anticipated  that  the 
solutions  will  be  mesh-size  dependent.  Calculations  are  run  for  different  crystal  sizes  and  several 
nucleation  densities.  In  all  instances,  the  simulations  are  simple  shear  with  the  single  slip  system 
orthogonal  to  the  direction  of  shear.  Periodic  boundary  conditions  are  applied  as  they  were  in 
section  3.1,  and  dislocation  transmission  is  prohibited  on  the  upper  and  lower  surfaces  where  the 
shear  is  applied.  The  yield  strength  is  50  MPa,  the  shear  modulus  30  GPa,  the  bulk  modulus 
60  GPa,  and  the  Burgers  vector  0.25xl0-3  pm. 
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4.2  Semi-discrete  Simulation  Results 


4.2.1  Low  Density  of  Nucleation  Sites 

Simple  shear  simulations  of  well-annealed  crystals  with  a  low  density  of  dislocation  nuclei  were 
run  for  model  regions  50  pm  wide  and  three  model  heights.  The  mesh  was  400  elements  wide. 
The  region  heights  were  5,10,  and  100  pm,  and  the  numbers  of  elements  in  the  thickness 
direction  were,  respectively,  40,  80,  and  800.  This  gives  an  element  size  of  0.125  pm  is  both  x 
and  y  directions  for  all  three  models.  The  number  of  dislocation  nucleation  sites  corresponding  to 
these  meshes  is  7,  14,  and  158,  respectively.  This  provides  roughly  the  same  number  of 

/T  _ r\ 

dislocations  nucleation  sites  per  unit  area  with  a  site  density  of  approximately  3.2x10  cm  .  This 
is  in  the  range  of  the  dislocation  density  for  a  well-annealed  metal. 

The  number  of  dislocations  that  passed  through  the  elements  at  a  shear  strain  of  0.01  is  presented 
in  figure  11.  The  complete  simulation  regions  are  shown  for  the  5-  and  10-pm-high  models.  Only 
the  top  and  bottom  10  pm  of  the  100-pm  model  are  shown,  since  most  of  the  center  section 
appears  as  lines  connecting  the  upper  and  lower  portions.  The  most  obvious  feature  is  the 
discrete  deformation.  Slip  occurs  only  along  slip  systems  containing  the  nucleation  sites.  The  slip 
traverses  the  crystal  vertically  along  lines  of  elements  that  contains  the  slip  planes.  The 
remaining  elements  are  elastic.  Each  of  the  7  nucleation  sites  produced  slip  in  the  5-pm  crystal, 
and  13  of  the  14  nucleation  sites  in  the  10-pm  model  are  evident.  One  of  the  nucleation  sites  in 
this  crystal  is  near  the  bottom  boundary  and  also  near  another  nucleation  site,  and  it  does  not  slip. 
The  100-pm  crystal  contains  more  nucleation  sites  and  many  active  slip  planes  are  evident. 
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Figure  1 1 .  Number  of  dislocations  passing  through  each  element  for  the  discrete  dislocation  simulations  of  three 
single  crystal  sizes  indicated.  The  color  levels  for  the  5-  and  10-pm  crystals  are  the  same.  The  center 
80  pm  is  omitted  from  the  100-pm  crystal  to  highlight  the  gradients  at  the  top  and  bottom  surfaces. 

An  important  feature  of  figure  1 1  is  the  slip  gradient.  Since  dislocations  cannot  pass  through  the 
upper  and  lower  boundaries,  the  slip  (number  of  dislocations  passed)  at  these  surfaces  is  zero. 
The  greatest  number  of  dislocations  has  passed  near  the  center  of  the  crystal.  The  slip 
distributions  along  the  13  active  slip  planes  of  the  10-pm  crystal  are  plotted  in  figure  12.  The 
average  of  the  13  curves  is  plotted  as  a  heavy  black  line,  and  the  individual  slip  results  are  shown 
as  the  thin  gray  lines.  All  of  the  curves  have  a  parabolic  appearance  near  the  boundaries.  Several 
have  a  flat  profile  across  the  center.  These  are  associated  with  slip  planes  that  have  dislocation 
nucleation  sites  near  the  surface,  and  they  are  also  near  more  highly  slipped  systems.  One  of  the 
dislocations  emitted  from  the  dipoles  travels  only  a  short  distance  to  the  boundary,  and  the  back 
stress  due  to  the  pile  up  reduces  the  driving  force  for  further  nucleation.  The  other  dislocation 
travels  much  of  the  way  across  the  crystal  before  it  meets  elevated  shear  resistance  from  the  pile- 
up  at  the  other  side  of  the  crystal.  The  long  travel  distance  creates  the  flat  center  region  of  the 
curve.  Also,  the  nucleation  sites  near  the  boundaries  produce  fewer  dislocations. 
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Figure  12.  Number  of  dislocations  passed  through  each  element  along  the  13 

active  slip  planes  in  figure  11  for  the  10-pm-high  crystal  plotted  against 
the  through  thickness  location.  The  light  lines  are  the  individual  results 
and  the  heavy,  dark  line  is  the  average. 

All  of  the  dislocations  must  lie  between  the  nucleation  sites  where  they  originate  and  the  crystal 
boundaries.  The  slip  distribution  is  directly  related  to  the  current  positions  of  the  dislocations. 

The  dislocation  positions  are  shown  in  figure  13.  As  with  figure  11,  the  entire  simulation  regions 
are  presented  for  the  5-  and  10-pm-thick  crystals,  and  only  the  upper  and  lower  10  pm  are  shown 
for  the  100-pm  thick  crystal.  Dislocations  of  opposite  sign  originating  at  the  dipoles  are  shown 
by  the  red  and  blue.  Dislocations  of  one  sign  move  to  the  top  and  those  of  the  other  sign  move  to 
the  bottom.  The  dislocation  density  is  greatest  at  the  boundaries  and  tapers  off  toward  the  center 
of  the  crystal.  This  is  the  classic  picture  of  an  edge  dislocation  pile-up.  The  dark  red  and  blue 
elements  on  the  crystal  interior  indicate  the  location  of  the  dipole  nucleation  sites  where 
dislocations  can  accumulate  before  gliding  toward  the  boundaries. 
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Figure  13.  Number  of  dislocations  currently  within  each  element  for  the  discrete  dislocation  simulations  of  three 
single  crystal  sizes.  The  color  levels  for  the  5-  and  10-pm  crystals  are  the  same.  The  center  80  pm  is 
omitted  from  the  100-pm  crystal  to  highlight  the  gradients  at  the  top  and  bottom  surfaces. 

Each  element  can  contain  multiple  dislocations,  and  the  elements  adjacent  to  the  boundaries  of 
the  smaller  two  crystal  sizes  contain  an  average  of  16.  That  average  drops  to  seven  in  the  next 
element  away  from  the  boundary,  demonstrating  a  large  gradient  in  dislocation  density  at  the 
crystal  boundaries.  The  first  six  elements  along  the  slip  planes  and  near  the  surface  all  contain 
dislocations  in  these  two  simulations.  Most  slip  planes  only  have  a  few  elements  without 
dislocations  along  their  entire  lengths.  Hence,  the  dislocation  density  is  relatively  dense  for  the 
5-  and  10- pm  crystal  heights. 

The  situation  is  different  for  the  100-pm-high  crystal.  With  many  more  nucleation  sites  and 
dislocations  to  carry  the  deformation,  the  driving  force  and  pile-ups  are  much  less  severe  for  the 
100-pm  crystal,  and  the  maximum  number  of  dislocations  per  element  is  four.  Here  the  gradient 
is  evident  in  the  sparseness  of  elements  containing  dislocations.  Near  the  surface  most  elements 
along  the  slip  planes  contain  dislocations,  but  as  the  distance  from  the  boundary  increases,  the 
elements  containing  dislocations  become  increasingly  sparse. 

The  discrete  slip  creates  a  highly  nonuniform  stress  field,  which  is  illustrated  by  the  in-plane 
stress  components  of  the  10-pm-high  crystal  presented  in  figure  14.  The  shear  stress,  Sig-XY,  is 
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low  along  the  slip  planes  due  to  sliding,  and  it  builds  up  between  the  slip  planes.  The  slip  planes 
effectively  partition  the  domain  into  finite  size  elastic  blocks,  each  undergoing  simple  shear  with 
the  shear  stress  on  the  lateral  boundaries  set  by  the  slip  resistance  on  the  slip  planes.  The  shear 
stress  rises  toward  the  interior  of  the  blocks  as  the  distance  from  the  shear  stress  lateral  boundary 
condition  increases.  Displacement  is  prescribed  on  the  upper  and  lower  boundaries.  This  imposes 
an  additional  constraint  that  keeps  the  boundary  nodes  uniformly  spaced  in  the  horizontal 
direction,  and  the  shear  stress  is  more  or  less  uniform  along  the  upper  and  lower  surfaces. 


Figure  14.  Plots  of  the  in-plane  stress  components  at  a  strain  of  0.01  for  the  10-pm-high 
single  crystal  deforming  by  discrete  single  slip  on  13  slip  planes. 

As  would  be  expected  from  an  elastically  sheared  block,  the  comers  experience  significant 
normal  stress,  Sig-YY,  as  the  blocks  try  to  shear  and  rotate  between  the  rigid  platens.  At  every 
material  point,  the  momentum  equation  in  the  vertical  direction  relates  the  gradient  of  Sig-YY  in 
the  vertical  direction  to  the  gradient  of  Sig-XY  in  the  horizontal  direction.  The  vertical  Sig-YY 
gradient  adjacent  to  the  slip  planes  is  small  near  the  boundary  and  becomes  approximately 
constant  (linear  stress  profile)  in  the  mid-thickness  of  the  crystal.  This  coincides  with  the 
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horizontal  shear  stress  gradient  adjacent  to  the  slip  planes  being  small  at  the  boundary  and  nearly 
constant  in  the  mid-thickness  region.  Likewise,  the  horizontal  gradient  in  the  Sig-XX  stress  is 
related  to  the  vertical  gradient  in  Sig-XY,  and  similar  gradient  correlations  can  be  made. 

4.2.2  Moderate  Density  of  Nucleation  Sites 

The  simulations  of  section  4.2.1  were  repeated  with  a  tenfold  increase  in  the  number  of 
dislocation  nucleation  sites.  These  were  92,  172,  and  1600  for  the  5-,  10-,  and  100-pm-high 
single  crystals.  This  gives  a  nucleation  site  density  of  approximately  3.2xl07  cm-2. 

Plots  of  the  dislocation  locations  are  shown  in  figure  15  for  the  5-  and  10-pm-thick  single 
crystals  deformed  in  simple  shear  to  a  strain  of  0.01.  As  in  the  simulation  with  a  low  nucleation 
site  density,  the  dislocation  dipoles  split  with  dislocations  of  one  sign  moving  toward  the  upper 
surface  and  dislocations  of  the  opposite  sign  move  toward  the  bottom  surface.  However,  since 
there  are  significantly  more  dislocations  to  accommodate  the  strain,  they  do  not  travel  as  far  as  in 
the  prior  analysis.  The  opposite  signed  dislocations  are  not  as  segregated  and  the  total  number  of 
dislocations  in  the  elements  at  the  boundaries  is  significantly  less.  The  results  for  the  100-pm 
simulation  are  not  shown.  This  calculation  had  even  more  nucleation  sites,  less  slip,  and  few 
dislocations  accumulated  at  the  boundaries.  Consequently,  segregation  of  the  opposite-signed 
dislocations  to  the  top  and  bottom  surfaces  was  not  apparent  at  a  strain  of  0.01. 


Figure  15.  Number  of  dislocations  within  each  element  for  discrete  dislocation  simulations  of  the  5-  and  10-|am 
single  crystals  with  a  moderate  density  of  nucleation  sites. 

The  dislocations  in  these  moderate  source  density  calculations  are  much  closer  together  than  in 
the  low  nucleation  source  runs  of  section  4.2.1.  As  a  result,  the  dislocations  are  too  close  for  the 
finite  element  discretization  to  resolve  the  stress  field  satisfactorily  (figure  16).  So  while  the 
motion  of  the  dislocations  in  figure  15  appears  to  be  rational,  the  lack  of  sufficient  spatial 
resolution  to  resolve  the  stress  gradients  renders  the  solution  questionable.  A  larger  element  size 
would  result  in  more  elements  containing  dislocations,  and  the  stress  resolution  would  degrade 
further.  Use  of  a  smaller  element  size  would  allow  more  elements  between  dislocations  to  better 
resolve  the  field  gradients. 
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Figure  16.  Shear  stress  distribution  for  lO-jam-thick  crystal  simulation  with  a  moderate 
dislocation  nucleation  site  density. 

4.2.3  Reduced  Element  Size 

Another  set  of  simulations  was  run  to  explore  the  effect  of  element  size.  The  mesh  dimensions 
used  in  sections  4.2.1  and  4.2.2  were  reduced  by  a  factor  of  5  in  each  direction  such  that  the 
crystal  width  is  10  pm  and  the  heights  were  1,  2,  and  20  pm.  The  element  dimension  for  these 
simulations  was  0.025  pm.  The  dislocation  nucleation  site  seeding  used  in  the  mesh  of  section 
4.2.2  was  also  used  here,  so  the  number  of  dislocation  nucleation  sites  is  the  same  as  in  section 
4.2.2.  However,  the  resulting  nucleation  site  density  is  25  times  greater  than  in  section  4.2.2 
(approximately  8x10s  cm-2  for  each  of  the  three  crystal  sizes)  because  the  mesh  is  a  factor  of  5 
smaller  in  each  direction.  An  alternative  configuration  would  have  been  to  keep  the  nucleation 
site  density  constant,  but  with  the  higher  spatial  resolution,  that  would  have  resulted  in  fewer  slip 
planes  than  are  evident  in  figure  11. 

The  dislocation  locations  from  these  three  simple  shear  runs  are  shown  in  figure  17  at  a  strain  of 
0.01.  As  with  the  previous  simulations,  the  full  crystal  heights  are  shown  for  the  smaller  two 
crystal  sizes  and  only  the  top  and  bottom  10%  of  the  tallest  crystal  is  shown.  Also  similar  to  the 
previous  simulations,  the  dipole  nucleation  sites  create  oppositely  signed  pairs  of  dislocations 
that  glide  up  and  down  toward  the  model  boundaries.  Pile-ups  are  evident,  and  it  appears  that  the 
dislocations  are  more  distinct  and  separated  relative  to  the  crystal  height  than  those  in  figures  13 
and  15.  This  is  expected  since  the  dislocation  spacing  is  a  physical  distance  and  the  model  has 
effectively  zoomed  in  by  a  factor  of  5. 
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Figure  17.  Number  of  dislocations  currently  within  each  element  for  the  discrete  dislocation  simulations  of  three 
single  crystal  sizes.  The  color  levels  for  the  1-  and  2-jam  crystals  are  the  same.  The  center  16  |am  is 
omitted  from  the  20-|am  crystal  to  highlight  the  gradients  at  the  top  and  bottom  surfaces. 


The  slip  planes  in  the  1-  and  2-pm-high  crystals  appear  distinct,  and  the  dislocation  pile-ups  lie 
cleanly  along  the  planes.  In  contrast,  several  of  the  nucleation  sites  in  the  20-pm-high  crystal 
have  horizontal  locations  similar  to  other  nucleation  sites,  resulting  in  closely  spaced  active  slip 
planes.  For  these  closely  spaced  slip  planes,  the  locations  of  the  dislocations  near  the  surface 
appear  to  be  coordinated.  Such  coordination  could  be  related  to  a  lower  energy  location  in  the 
interacting  stress  fields. 

The  shear  stress  field  for  this  higher  resolution  simulation  is  shown  in  figure  18.  As  with  the 
shear  stress  shown  in  figure  16,  the  spatial  discretization  is  too  coarse  compared  to  the 
dislocation  spacing  to  resolve  the  gradients.  Many  of  the  lowest  stresses,  indicated  by  blue  in 
figure  18,  are  isolated.  Consequently,  the  quality  of  the  solutions  is  questionable. 
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Figure  18.  Shear  stress  distribution  for  the  10  |am  x  2  |am  crystal  in  a  calculation  with  an 

element  size  of  0.025  |jm  and  172  dislocation  nucleation  sites.  The  stress  gradients 
are  under-resolved. 


An  important  piece  of  information  on  the  shear  stress  plots  is  the  minimum  and  maximum  stress. 
In  figure  16,  the  maximum  is  120  MPa  and  the  minimum  is  -18  Mpa,  whereas  these  are, 
respectively,  604  MPa  and  -307  MPa  in  figure  18.  The  reason  can  be  traced  to  the  quantized 
shear  computed  through  equations  12  and  13.  As  a  dislocation  is  about  to  leave  an  element,  it 
exerts  an  elastic  shear  strain/stress  on  the  neighboring  element.  As  the  dislocation  crosses  the 
element  boundary  and  traverses  the  element,  the  shear  stress  through  goes  to  zero  and  reverses 
sign.  The  stress  due  to  the  discrete  shear  peaks,  and  it  has  the  opposite  sign  as  the  dislocation  is 
leaving  the  element.  Analytically,  the  shear  stress  is  infinite  at  the  dislocation  core,  and  the 
spatial  discretization  effectively  averages  the  stress  over  the  element  size.  As  the  element  size 
gets  smaller,  the  peak  stress  represented  in  the  mesh  will  increase,  but  the  gradient  is  less  well 
resolved.  In  figure  18,  the  peak  shear  stresses  are  an  order  of  magnitude  higher  than  the  stress 
required  for  slip.  Such  high  stresses  dominate  the  solution  and  the  error  in  the  gradient 
overwhelms  the  slip  resistance  of  the  crystal.  This  is  another  reason  that  the  quality  of  the 
solution  is  questionable. 

4.2.4  Stress-strain  Response  from  Semi-discrete  Model 

The  shear  stress-strain  curves  from  the  semi-discrete  dislocation  simulations  are  shown  in 
figure  19.  The  element  size,  nucleation  site  density,  crystal  width,  and  crystal  aspect  ratio  are 
indicated  in  the  legends.  The  results  from  all  of  the  simulations  in  section  4.2  are  presented  in 
addition  to  a  set  of  three  runs  similar  to  those  of  section  4.2.2  but  with  the  domains  and  elements 
10  times  larger.  As  discussed  in  earlier,  the  dislocations  in  this  calculation  are  significantly 
under-resolved. 
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Shear  Strain 


Figure  19.  Shear  stress-strain  response  for  the  semi-discrete  crystal  simulations  of  section  4.  The  legends  and  arrows 
show  the  model  width,  crystal  aspect  ratio,  element  length,  and  the  percent  of  elements  containing 
dislocation  nucleation  sites  along  with  the  corresponding  nucleation  site  density.  Curves  of  the  same 
color  have  the  same  element  size  and  element  site  nucleation  percentage  but  different  model  aspect 
ratios. 


The  color  groupings  contain  sets  of  three  simulations  where  only  the  crystal  height  was  varied. 
Within  these  sets  of  three  curves  of  the  same  color,  the  results  from  the  short  crystals  are 
represented  by  dotted  lines,  the  intermediate  height  crystals  are  indicated  by  the  dashed  lines,  and 
the  tall  crystals  by  the  solid  lines.  In  all  cases,  the  stress  increases  as  the  crystal  height  decreases. 

The  green  curves  show  the  results  from  section  4.2.1  where  the  dislocation  nucleation  sites  were 
sparse.  The  short  crystal  contained  only  seven  active  slip  planes,  and  relatively  large  elastic 
regions  are  providing  the  shear  stress  resistance.  The  few  slip  planes  do  little  to  relieve  the 
stress,  and  the  result  is  nearly  elastic.  The  intermediate  height  crystal  was  analyzed  in  some 
detail  above.  The  additional  slip  planes  created  smaller  elastic  domains,  and  the  stress  relief 
provided  by  the  slip  planes  is  sufficient  to  have  an  appreciable  impact  on  the  overall  stress.  With 
the  additional  nucleation  sites  in  the  tall  crystal,  the  elastic  domains  are  small  and  the  shear  stress 
peaks  and  then  saturates  after  a  small  drop. 
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The  nucleation  site  density  was  increased  by  an  order  of  magnitude  over  the  green  results  to 
produce  the  red  results,  section  4.2.2.  The  mesh  and  model  size  are  identical.  Here  there  are 
sufficient  dislocations  to  allow  pervasive  slip  even  for  the  shortest  crystal.  The  stress  saturates 
and  is  nearly  the  same  for  all  three  crystal  heights.  Using  the  same  number  of  dislocation  sites, 
the  mesh  size  is  then  increased  by  a  factor  of  10  in  both  directions  over  the  red  curves  to  produce 
the  black  curves.  Here  too,  the  stress  saturates  and  simulations  for  all  three  crystal  heights  give 
approximately  the  same  shear  stress.  The  shear  stress  is  lower  than  the  red  curve  because  the 
singular  shear  stress  near  the  dislocations  is  averaged  over  a  large  volume,  and  the  resulting 
lower  element  stress  near  the  dislocation  lowers  the  overall  shear  stress.  This  simulation  was  not 
described  in  earlier  sections. 

The  calculations  producing  the  blue  curves  are  described  in  section  4.2.3.  This  simulation  had  a 
mesh  a  factor  of  5  smaller  in  each  direction  than  the  red  curve,  and  the  number  of  dislocation 
nucleation  sites  was  the  same.  Even  though  there  appears  to  be  a  sufficient  number  of 
dislocations  in  figure  17  to  produce  a  smooth  plastic  response,  the  shear  stress-strain  curves  are 
very  noisy.  The  reason  is  the  very  high  stress  near  the  dislocations  that  results  from  the  finer 
spatial  discretization.  Even  though  stress  jumps  occur  over  small  volumes  as  the  dislocations 
move  from  one  element  to  another,  the  magnitude  of  the  jump  is  sufficient  to  affect  the  average 
stress. 


5.  Discussion 


The  results  from  the  continuum  slip  gradient  model  presented  in  section  3  demonstrate  that  an 
operator  split  approach  can  be  used  with  an  explicit  dynamic  time  integration  scheme  to  include 
gradient  effects  in  crystal  plasticity  simulations.  The  nonlocal  gradient  formulation  was  an 
outgrowth  of  a  penalty  approach  to  enforce  slip  continuity  between  neighboring  elements.  The 
model  was  implemented  into  a  large-scale  parallel  code  with  minimal  disruption  to  the  flow  of 
the  calculations.  There  appear  to  be  additional  time  step  restrictions  for  smaller  element  sizes, 
but  the  nature  of  the  restrictions  was  not  explored  in  the  current  investigation. 

The  single  crystal  and  polycrystal  results  show  the  anticipated  diffusion  of  sharp  deformation 
fields  and  the  expected  stress-strain  trends  with  size  scale.  While  having  the  appearance  of  a 
successful  modeling  effort,  a  critical  comparison  of  the  model  assumptions  to  the  physical 
configuration,  which  it  is  supposed  to  represent,  reveals  a  significant  disconnect.  This  is 
particularly  evident  in  the  micron  size  range  as  the  spatial  resolution  of  the  finite  element  grid  is 
refined.  An  inherent  assumption  in  the  finite  element  crystal  plasticity  model  is  that 
microstructure  features  are  either  explicitly  resolved  by  the  grid  or  the  microstructure  is  at  a  scale 
sufficiently  smaller  than  the  element  size  such  that  it  can  be  approximated  as  a  smoothly  varying 
field.  Similarly,  a  nonlocal  gradient  extension  to  the  crystal  model  should  resolve  the  gradient 
over  several  elements.  As  indicated  above  in  table  1,  the  element  size  must  be  larger  than  1  pm 
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to  have  an  average  of  one  dislocation  per  element  for  a  well-annealed  material.  At  lower 
resolutions,  the  field  within  the  element  is  not  smooth  and  the  discreteness  of  the  microstructure 
is  not  resolved.  A  continuum  crystal  model  at  this  or  finer  resolution  does  not  provide  a  valid 
representation  of  the  physical  microstructure  for  an  annealed  crystal.  Consequently,  simulations 
cannot  be  expected  to  capture  the  microstructure  mechanisms  at  this  scale  accurately. 

Looking  strictly  at  the  dislocation  density  for  heavily  deformed  metals,  from  table  1,  it  appears 
that  the  crystal  model  with  submicron  spatial  resolution  could  represent  a  smeared  dislocation 
microstructure  within  each  element.  However,  dislocations  are  rarely  uniformly  distributed  in 
heavily  deformed  metals.  They  typically  organize  into  walls,  which  create  a  cell  structure.  The 
cell  walls  have  a  very  high  dislocation  density,  and  the  cell  interiors  have  a  low  dislocation 
density.  Extensive  analysis  of  heavily  deformed  nickel  by  Hughes  and  Hansen  (2000)  shows  that 
these  cell  sizes  are  greater  than  0.1  pm  after  cold  rolling  to  a  98%  reduction,  or  a  logarithmic 
strain  in  excess  of  3.9.  In  order  to  have  a  smeared  representation  of  such  a  cell  microstructure 
within  each  element,  the  element  size  would  have  to  be  on  the  order  of  1  pm  or  larger. 
Traditional  crystal  plasticity  simulations  using  smaller  element  sizes  will  disregard  a 
microstructure  that  could  be  resolved  at  that  spatial  resolution,  and  the  computed  deformation 
mechanisms  may  not  be  accurate. 

As  in  the  modeling  of  composites,  refinement  of  the  spatial  discretization  for  the  standard 
continuum  crystal  plasticity  models  has  a  limit  when  the  element  size  begins  to  approach  the 
length  scale  of  important  microstructure  features.  In  order  to  have  predictive  microstructure- 
based  simulations,  the  modeling  strategy  has  to  change  with  mesh  refinement  to  match  the 
physical  microstructure  at  that  refinement  level.  When  the  element  size  during  mesh  refinement 
of  a  composite  approaches  the  size  of  the  reinforcement,  a  homogenized  continuum  model  is  no 
longer  appropriate,  and  the  reinforcement  and  matrix  phases  should  be  modeled  explicitly.  There 
is  nothing  to  prevent  analysis  using  a  homogenized  model  with  elements  smaller  than  the 
reinforcement,  but  the  real  mechanisms  operating  at  that  scale  would  be  missed  entirely.  The 
calculations  would  not  be  predictive  of  mechanisms  occurring  at  that  resolution. 

For  pure  metals,  the  limiting  microstructure  features  for  a  crystal  plasticity  model  are  the 
dislocation  spacing  or  the  dislocation  cell  size,  whichever  is  larger.  These  are  typically  on  the 
order  of  0.1  to  1.0  pm.  Further  mesh  refinement  near  this  range  of  element  sizes  requires  some 
explicit  accounting  for  the  dislocation  microstructure  and  a  change  in  the  underlying  constitutive 
relations.  The  continuum  crystal  model,  including  gradients,  will  run  with  smaller  elements,  but 
the  solution  will  miss  the  effects  of  the  discrete  dislocation  interactions  and  any  mechanisms 
specific  to  dense  dislocation  walls  and  cell  structures.  It  would  give  a  false  sense  of  a  high 
resolution  simulation  since  the  microstructure  governing  the  behavior  is  not  captured  spatially. 
The  importance  of  accounting  for  the  discreteness  of  dislocations  in  a  pile  up  is  highlighted  by 
Roy  et  al.  (2008)  and  Baskaran,  et  al.  (2010).  In  general,  however,  mechanisms  at  this  length 
scale  have  received  little  attention.  Modeling  at  the  size  scale  of  precipitate  or  phase 
microstructure  in  alloys  would  also  require  further  length  scale  considerations. 
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Nonlocal  gradient  calculations  implicitly  rely  on  the  validity  of  the  representation  of  the 
microstructure  within  the  elements.  If  the  element  size  is  of  the  order  of  important  microstructure 
features  such  that  deformation  mechanisms  are  not  represented  appropriately,  the  gradient 
calculations  would  not  be  based  on  reliable  element  data.  Hence,  gradient  calculations  are 
restricted  to  the  same  lower  resolution  limits. 

A  semi-discrete  model  based  on  the  standard  crystal  plasticity  model  was  investigated  as  an 
expedient  method  of  crudely  accounting  for  the  discrete  nature  of  dislocations  at  finer  spatial 
resolutions.  Only  elements  containing  dislocations  could  slip  plastically,  and  dislocations 
produce  a  definite  level  of  slip  after  passing  through  an  element.  While  the  results  of  this  simple 
model  showed  dislocation  pile-ups  and  other  expected  features,  the  specific  approach  was 
unsatisfactory.  The  overriding  issue  is  the  singular  stress  field  from  the  dislocations  dominating 
the  solution  as  the  mesh  was  refined  to  resolve  the  stress  gradients.  The  stress  field  is  driven  by 
the  quantization  of  slip,  which  provides  a  higher  stress  magnitude  related  to  the  better  resolution 
of  the  singularity  at  finer  spatial  resolutions.  There  does  not  appear  to  be  a  range  of  element  sizes 
that  will  both  resolve  the  stress  field  and  not  suffer  from  the  effects  of  the  singularity.  Perhaps 
other  semi-discrete  approaches  could  be  successful;  this  one  was  not. 

Discrete  dislocation  dynamics  simulations  (e.g.,  Kubin  et  al.,  1992  and  Arsenlis  et  al.,  2007) 
explicitly  account  for  dislocations  and  their  interactions  and  provide  one  means  for  incorporating 
dislocation  microstructure  at  finer  spatial  resolutions.  Finite  element  methods  have  been  coupled 
with  the  discrete  dislocation  simulations  though  several  approaches  (Van  der  Giessen  and 
Needleman,  1995;  Fivel  et  al.,  1998;  Yasin  et  al.,  2001).  These  types  of  formulations  may  be 
employed  at  the  finer  spatial  resolutions  in  multiscale  modeling  schemes.  There  could  be  a 
transition  from  traditional  continuum  crystal  plasticity  to  such  a  representation  when  the  spatial 
resolution  is  fine  enough,  as  in  Wallin  et  al.  (2008). 

An  area  that  has  seen  little  model  development  activity  is  deformation  associated  with 
dislocation  cell  walls.  In  addition  to  being  sources  and  sinks  of  dislocations  for  the  cells, 
dislocations  may  run  within  some  cell  walls,  causing  slip.  These  walls  could  be  treated  as 
entities  in  resolved  calculations  of  dislocations  cells  where  they  would  interact  with  models  of 
the  cell  interior  which  would  transmit  dislocations  across  the  cells.  This  approach  may  be  more 
efficient  than  modeling  every  dislocation  within  the  cell  walls.  This  is  an  area  for  future  research. 


6.  Conclusion 


A  nonlocal  crystal  model  based  on  the  second  gradient  of  crystal  slip  was  implemented  in  a 
large-scale  parallel  finite  element  code,  and  the  results  show  the  expected  trends  of  decreasing 
the  severity  of  gradients  and  increasing  strength  with  decreasing  physical  size.  Although  the 
model  appears  successful  on  the  surface,  there  are  concerns  over  the  adequacy  of  the  traditional 
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crystal  model  to  represent  the  microstructure  at  smaller  size  scales.  Elements  smaller  than  the 
dislocation  spacing  or  the  dislocation  cell  dimensions  cannot  capture  features  of  the  deformation 
at  that  size  scale  accurately.  This  prompted  exploration  of  a  semi-discrete  crystal  plasticity 
model.  The  semi-discrete  model  also  showed  realistic  results  in  terms  of  dislocation  pile-ups,  but 
the  stress  solution  at  the  dislocation  singularity  dominated  the  solution  at  smaller  element  sizes, 
rendering  the  solution  noisy  and  very  mesh  dependent. 

The  size  scale  of  the  microstructure  of  a  typical  pure  metal,  in  terms  of  dislocation  spacing  and 
size  of  dislocation  cells,  is  on  the  order  of  0.1  to  1.0  pm.  If  the  element  size  from  a  simulation  is 
small  enough  to  be  within  this  range,  the  microstructure  features  and  mechanism  should  be 
modeled  explicitly  in  order  to  capture  the  lower  length  scale  behavior.  Models  particular  to  this 
size  scale  are  needed. 
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